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CALCULATION OF BOUNDARY LAYERS OF OSCILLATING AIRFOILS 
Tuncer Cebeci* and Lawrence W. Carr** 

ABSTRACT 

A two-point finite difference unsteady laminar and turbulent 
boundary layer computational method has been used to investigate 
the properties of the flow around an airfoil (NACA 0012) oscillating 
through angles of attack up to 18 degrees, for reduced frequencies 
of 0.01 and 0.20. The unsteady potential flow was determined using 
the unsteady potential flow method of Geissler. The influence of 
transition location on stall behavior was investigated, using both 
experimentally determined transition information, and transition 
located at the pressure peak; the results show the need for viscous- 
inviscid interaction in future computation of such flows. 


*Mechamcal Engineering Dept., California State University, Long 
Beach, CA 90840. 

**Ames Research Center and Aeromechanics Laboratory, AVRADC0M 
Research and Technology Laboratories, Moffett Field, California. 



INTRODUCTION 


Present knowledge of steady boundary-layer flows is considerable and 
stems from exhaustive investigation over almost eighty years. In contrast, 
unsteady boundary-layer flows have received slight attention and only in the 
recent past. The requirements imposed by the need to improve the performance 
of helicopter rotors, wind-energy devices and hydrofoil characteristics imply 
tne need for further research. This is particularly so since many of the 
unsteady flow investigations have been concerned with detail such as the 
existence or the nonexistence of singularity and its structure, see for example 
the papers contained in ref. 1. Here we are concerned with aspects of our 
ability to calculate the flow properties around oscillating airfoils in a range 
of parameters of direct relevance to engineering practice. 

The need for further development of boundary -layer procedures to represent 

9 

unsteady flows is clear from the recent experiments of McCroskey et al. 1 ', 

Tijdeman~ l , Davis and Malcolm 4 , Cousteix et al. 5 , Carr et al. 5 . Young 7 , Carr 
8 9 

and McAlister and Geissler . The range of measurements is reasonably extensive 
and includes subsonic and transonic flow, a range of airfoil sections and 
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parameters of unsteady motion. Four flow regimes have been identified^ and 
correspond to no stall, stall onset, light stall and deep stall. None of 
these regimes nor the apparent breakdown of the unsteady boundary layer to 
form a large vortex near the surface at large angle of attack have been 
represented satisfactorily by calculation methods. 

Previous attempts to calculate oscillating airfoil flows have involved 

11 12 

either the solution of Navier-Stokes equations, Mehta and Shamroth or 

13 14 

the solution of the boundary-layer equations, Cebeci and Carr * and 
1 5 

Geissler . While both approaches have merits, it is clear that the latter 
procedure, with further development, is likely to allow more precise solution 
of the equations and economy of computer resource:-. 

The calculation of boundary-layers on an oscillating airfoil differs from 
the usual nonoscillating airfoil-flow calculations in that difficulties are 
caused by the translation of the stagnation point in space and time. In par- 
ticular, it is necessary to develop a procedure to generate initial conditions 
in the immediate vicinity of the moving stagnation point and to account for 
the flow reversal that occurs in this region. In reference 14, the present 
authors described two procedures for generating initial conditions near the 
stagnation point of an oscillating airfoil. The first procedure used the 
characteristic box scheme and was shown to be accurate and free from the limi- 
tations which may be imposed by the flow reversal . The second procedure used 
a quasi-steady approach and was also shown to be appropriate provided that the 
region where it was used was far from any flow reversal . 

The work described here is the continuation of that of reference 14, and 
explores the calculation of laminar and turbulent boundary layers on the whole 
oscillating airfoil. It has three separate but related phases. First, it is 
necessary to conduct numerical tests of the procedures of ref. 14 and of the 
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boundary conditions generated by the potential flow code of ref. 16 . 

Secondly, oscillating airfoils involve transition from laminar to turbulent 
flow and the data provided by Carr and his coworkers allow the direct 
evaluation of transition assumptions which must be made in the boundary-layer 
calculations. In particular, a simple assumption is to involve transition at 
the location corresponding to the peak pressure calculated by the inviscid 
flow method. The available experimental data permit the implications of this 
assumption to be evaluated and compared with those empirical transition 
methods used for steady flows and described in ref. 17. Thirdly, leading edge 
separation bubbles and trailing edge separation can also occur at higher 
angles of attack as demonstrated by Carr et al. S The ability of a 
calculation method to represent these regions of separation and the 
expressions used to model turbulent flows at low Reynolds numbers with 
separation as well as transitional flows need to be investigated. 

The study reported here is continuing and this contribution may be 
regarded as a report of progress since that of ref 14. 


BASIC EQUATIONS 


Boundary-Layer Equations 

The boundary-layer equations for an incompressible laminar or turbulent 
flow on an oscillating airfoil are well known and, with the eddy viscosity 
concept, can be written as 
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Here x denotes distance along the surface of the airfoil, y along the 
normal and b = v + e m . In the absence of mass transfer, Eqs. (1) and (2) are 
subject to boundary conditions given by 

y = 0; u = v = 0: y = <5, u = u e (x,t) (3) 

The presence of the eddy viscosity e m requires a turbulence model and we 
use the algebraic eddy-viscosity formulation developed by Cebeci and Smith 
(ref. 18). According to this formulation is defined by two separate 
formulas. In the inner region of the boundary layer (e m )-j is defined as: 


(e m ) i = {0.4y[l - exp{-y/A)]} 2 ||y-| y tr 0 < y 1 y c 
where 
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In Eq. (4), y^ r is an i ntermi ttency factor t'.iat accounts for the transi- 
tional region that exists between a laminar and turbulent flow. It is defined 
by 


Y tr 


- exp [-G (x - x t ) 




( 6 ) 


Here x t is the location of the start of transition and the empirical factor 
G which has the dimensions of velocity/( length) ', is given by (ref. 18) 


r 1 u e D -1 .34 

b " TZOU T K x. 

v tr 
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The transition Reynolds number is defined as R = (u x/v) t . 

x tr e 

In the outer region (e m )o 1s defined by: 
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The boundary between the inner and outer regions, y , is established by the 

V 

continuity of the eddy-viscosity formulas. 


Initial Conditions 

If initial conditions in the (t,y) plane are given at a station x Q on the 
upper surface of the airfoil and satisfy the condition u > 0 and, in addition, 
initial conditions are given in the (x,y) plane at t = 0, then the solution of 
Eqs. (1), (2) and (3) may be obtained for x > x Q and t > 0 until they 
breakdown (flow separation). A similar remark applies to the lower surface 
except that u < 0. The initial conditions at t = 0 can be generated for both 
surfaces if steady conditions are assumed to prevail at that time. It is only 
necessary to solve the appropriate equations which, in this case, are given by 
Eq. (1 ) and by 


u vM= u +1. (b^) 
ax ay e djF ay v ay 7 


(9) 


There is no problem with the initial conditions for Eqs. (1) and (9) since the 
calculations start at the stagnation point x = x $ , v/h ere u g and u are zero 
for all y. 

Unlike steady flows, where u e and u are zero for all y at the 
stagnation point, the stagnation point is not fixed in an unsteady flow; 
although u g is zero, we cannot assume a priori that u is also zero. Ue nay 
avoid these difficulties by using an implicit method, but now we are faced 
with the problem of generating a starting profile on the new time line. A 
convenient and accurate procedure to calculate the first velocity profile at 
the new time-line has been developed by the present authors as described in 
ref. 14; it involves the use of the Characteristic Box scheme developed by 
Cebeci and Stewartson (ref. 19). Another procedure is to use a quasi -steady 
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approach in the immediate vicinity of the stagnation point. As long as the 
region where this approach is used does not contain flow separation, it is 
simpler to use than the first approach which employs the characteristic-box 
scheme. As a result, it is used here. 


Transformed Equations 

As in previous studies (see, for example, ref. 17), we use similarity 
variables to transform the governing equations before we seek their solution. 
For a steady flow, we use the Falkner-Skan transformation defined by 


n - / 



i> = /u gV x f(x,n) 


( 10 ) 


where \|> is the usual definition of stream function that satisfies the conti- 
nuity equation, Eq. (1), that is 


= 3£ v _ ii 
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With this transformation, Eq. (9) and its boundary conditions, Eq. (3), can be 
written as 


(bf")' + ff" + n[l - (f 1 ) 2 ] = x(f ' §£- - f" ||) (12) 

n = 0; f = f' = 0: n = n e ; f' = l (13) 


where primes denote differentiation with respect to ri, m denotes a dimen- 
sionless pressure- gradient parameter defined by 


and b is defined by 
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For unsteady flows, we use a transformation similar to that defined by 
Eq. (10) except that u g is now a function of both x and t and the dimension- 
less stream function F is a function of x, t and ^ ; we let 


u p (x,t) 

£ = / - v - x - — y. = /u e (x,t)v'x F(x,t,^) (15) 


With this transformation it can be shown that the continuity and momentum 
equations and their boundary conditions for unsteady incompressible flows can 
be written as 
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Solution Procedure 

Ue use Keller's Box method to solve the governing equations of the previ- 
ous section. This is a two-point finite-difference method which has been used 
to solve a wide range of parabolic partial -differential equations as discussed 
in ref. 19. The solution procedure for Eqs. (12) and (13) is identical to 
that described in ref. 17. The solution procedure used to generate the 
initial conditions in the (x,y) and (t,y) planes is described in reference 
13. For unsteady flows, where we now solve Eqs. (16) and (17), we use the 
solution procedure described in ref. 20. In regions where there are no flow 
reversals across the layer, we use the Standard Box scheme and in regions 
where there is flow reversal, we use the Zig-Zag scheme. 
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RESULTS AND DISCUSSION 


To accomplish the objectives stated in the Introduction with the 
boundary-layer method discussed in the previous section, we have considered 
the NACA 0012 airfoil and calculated its pressure distribution with Geissler's 
inviscid code 16 for two reduced frequencies, k = 0.01 and k = 0.20 and for 
two angles of attack 

f 5° + 5° sirtit (19a) 

1 8° + 10° sinwt (19b) 

The experimental data taken for these two angles of attack and two frequencies 

show that the one with the smaller variation in angle of attack, Eq. (19a), 

has no trail ing-edge separation on the whole airfoil for either frequency. 

Its maximum angle of attack range falls in what McCroskey and Pucci call the "no 

stall region" ( a _.„ < 13°). The only flow separation occurs in the form of a 
max 

bubble near the upper surface leading edge. For a > 5°, flow separates and 
reattaches as the flow goes from laminar to turbulent. 

The flow corresponding to the larger variation in angle of attack, Fq. 
(19b), on the other hand, has a leading edge separation bubble as well as an 
open trailing edge separation. Its maximum angle of attack range falls in 
between light stall (ot m ,„ <15°) and deep stall (« „ < 20°) with the 

magnitude of reduced frequency playing an important role in the performance of 
airfoil characteristics. As discussed by McCroskey and Pucci 10 , the quali- 
tative behavior of light stall in this case is sensitive to reduced frequency 
and maximum incidence for a specified airfoil at zero Mach number. The 
qualitative behavior is closely related to the boundary-layer separation 
characteristics (leading-edge vs trail ing-edge separation) and to the changes 
in this separation behavior a max and k. An important point that they note 
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is the amount of trai ling-edge separation suppressed by increasing the 
frequency. As a result of this unsteady effect, an airfoil that falls in the 
light stall region can pass from light-stall with k £ 0.10 to stall onset 
(a =14) for k > 0.15. 

Itlo X 

The boundary-layer method described in the previous section is able to 
calculate both laminar and turbulent flows. The calculations start at the 
stagnation point for the given pressure distribution and for the specified 
freestream. The initial conditions at t = 0 are obtained by solving the 
steady flow form of the equations and those in the (t.y) plane with the 
quasi-steady approach described earlier and applied only in the immediate 
neighborhood of the stagnation region. Transition is achieved by specifying 
its location as part of the input to the computer program. In contrast to 
steady flows where empirical correlations for transition have been developed 
and compared satisfactorily with experimental data, little is known of the 
same phenomena in unsteady flows. Here we have taken two approaches. In the 
first, if available, transition location has been specified directly from 
experiment and in the second it has been specified at the location of maximum 
pressure. The first assumption allows the evaluation of the numerical 
technique and the second is expected to be physically sound since the pressure 
distribution becomes increasingly peaky as the angle of attack is increased 
and the merits of this assumption can be tested. It can also be expected that 
the second assumption will become less appropriate as the peak in the pressure 
distribution diminishes with decreasing angle of attack. 

Figures 1 and 2 show the variation of displacement thickness 5 *, and 
local skin-friction coefficient, c^, distributions on the airfoil at 
different angles of attack for two reduced frequencies. Here <s* and c^ 
are defined by 
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(20a) 
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and in terms of the transformation given by Eq. (15) can be expressed as 


5* --zz (n* - FJ. 
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w 
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Transition was assumed to correspond to maximum pressure and the chord 
Reynolds number, R , was taken as 3 x 10^ and approximately 100 stations 
were taken in the x-direction with 41 stations in the time direction 
corresponding to Au>t = 9°. The steady-state conditions were obtained from 
the pressure distribution that corresponded to a ■■ 0 by taking cot = 270° and 
the calculations for t > 0 were performed for increasing angle of attack. Since 
we are using an unsteady pressure distribution to calculate the initial conditions 
at t = 0 and au/3t is not zero at t = 0, slight oscillation occurs when 
unsteady-flow calculations are performed for the next time step. For this 
reason a smoothing procedure was applied to the calculated results at that 
time station and the next few stations. 

As can be seen from the results, there is no flow separation on the 
airfoil except for a very small region at the trailing edge. This is without 
doubt due to the inaccuracy of the inviscid pressure distribution in that 
region where the effect of viscous forces on the inviscid flow is greatest. 

An inviscid-viscous interaction should eliminate the trail ing-edge separation. 

The calculations for case 2 are more difficult to perform than those for 
case 1 since the flow regime falls between light stall and deep stall, due to 
larger variations in angle of attack. Whereas there is almost no flow 
separation in case 1, there is both leading-edge and trailing-edge separation 
for case 2. Once the calculations have started for steady state for, say cot 
corresponding to 270°, we expect that flow separation will take place as the 
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angle of attack increases. However, whether the first flow separation will 
manifest itself in the form of leading-edge separation or trailing-edge sepa- 
ration remains to be explored. In either case, if the flow separates at the 
leading-edge at higher angles of attack, then the use of the boundary-layer 
method will lead to breakdown of the solutions due to the singularity at sepa- 
ration, and, as a result, calculations cannot be performed for the downward 
stroke beyond the separation point. For the case when the separation takes 
place near the trailing edge, again the calculations cannot be extended beyond 
the separation point for the downward stroke. 

The calculations for case 2 were again performed for the same chord 
Reynolds number and for the same number of x- and t-stations. The initial 
calculations were started by using the pressure distribution which corresponded 
to a = -2°. Transition locations were specified by two different approaches 
for two reduced frequencies, k = 0.01 and k = 0.20. 

Figure 3 shows the variation of transition as a function of angle of 
attack for case 2 at two frequencies. We see from these figures that the 
effect of reduced frequency on the location of transition is considerable. At 
the lower frequency, transition occurs almost at the same location on the 
upward and downward strokes whereas, with the higher frequency, the transition 
locations on both strokes differ noticeably. Compared with the lower frequency 
case, with increasing frequency, the transition location moves backward in the 
upward stroke and forward in the downward stroke. It is obvious that a correct 
prediction of airfoil characteristics in the same angle of attack range for 
different frequencies will depend on our ability to predict transition as a 
function of frequency among other things. 

Figures 4 to 8 show the results for case 2. Those in Figures 4 and 5 
correspond to a transition location at the maximum pressure peak. The first 
traces of flow reversal occur near the trailing edge at a * 8° and are 
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limited to less than the chord length. As the angle of attack increases, the 
flow reversal region near the trailing edge slightly increases becoming 
approximately 5% at a = 18°. There is no leading-edge separation for either 
frequency. 

Figure 6 shows the results for case 2 in which transition locations were 
input from the experimental data, according to Fig. 3b. Again the calculations 
were started at the same angle of attack, a = -2, and continued in the upward 
stroke mode. In this case, however, the solutions broke down at a = 11.09°. 
There was essentially no flow reversal and no signs of numerical difficulties 
until a = 9.54 when flow separation appeared around the trailing edge. At 
the next angle of attack, a = 11.09, a laminar leading-edge separation bubble 
appeared at around 4% chord and the solutions broke down shortly thereafter. 
Figure 6 also shows that while the computed displacement thickness distribu- 
tions are smooth, those corresponding to local skin friction are not. The 
wiggles in the latter case are a result of transition locations specified at 
different angles of attack. For example, we see from Figure 3b that for 
-2 < a < 0, the transition location moves from (x/c) tr = 0.40 to approximately 
(x/c) tr = 0.55. This means that the laminar flow calculations for a = 0 
contain a combination of laminar and turbulent flow characteristics which 
originated at a = -2 where the laminar flow calculations terminated at 
(x/c) tr = 0.40. 

Calculations were also performed with the transition criterion of Michel, 
which was devised for steady flows, and the results were found to be similar 
to those of Figure 6 except that the solutions broke down earlier at a =8°. 
Figure 7 allows comparison between the measured transition locations and those 
calculated from Michel's formula during the upward stroke together with those 
corresponding to maximum pressure peak. It is clear, therefore, that the use 
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of this correlation, although imperfect, is to be preferred to the specifica- 
tion of transition at the location of the pressure peak. 

The intermittency expression used for the transitional region in the eddy- 
viscosity formula can play an important role in the calculation of leading-edge 
separation bubbles. It is likely that this formula, which has been used 
satisfactorily for attached flows undergoing transition from laminar to turbu- 
lent flows, may not be as accurate for flows undergoing transition due to 
separation. It is evident from Eqs. (6) and (7) that the coefficient of the 
parameter G controls the transition length. Thus a decrease in G by a factor 
of 10 will reduce the transition length which, in turn, will lessen the pos- 
sibility of separation. Calculations shown in Figure 8 performed with a value 
of 1/120 for the coefficient of G revealed that the leading-edge separation had 
been suppressed allowing the calculations to proceed and to reveal trai ling- 
edge separation. This procedure permitted the calculations to be performed in 
the upward stroke to an angle of 18° where the trai ling-edge separation had 
moved forward to around 20% chord. This sensitivity of the calculated flow to 
the coefficient of G requires further confirmation which may be obtained, in 
part, by the use of interactive methods; in this way the breakdown caused by 
the singularity will no longer occur. 


CONCLUDING REMARKS 


The numerical method of ref. 13 has been used to investigate the proper- 
ties of the flow around an airfoil (NACA 0012) oscillating with angles of 
attack up to 18° and for reduced frequencies of 0.01 and 0.20. This encompas- 
ses the regions of no-stall and deep stall identified by McCroskey and 
Pucci^. The inviscid pressure distribution was obtained from Geissler's 
method. The calculations were performed without numerical difficulties in 
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both flow regimes and the zig-zag numerical scheme was shown to work satis- 
factorily in regions of flow reversal. 

The influence of transition location was investigated with the numerical 
method. The specification of transition at the location of pressure peak led 
to results which were in disagreement with experiment for the low frequency 
case; in particular, the calculated results did not reveal the deep stall 
observed in ref. 10. In contrast, the results for the high frequency case were 
in accord with experiment. Use of the experimentally determined transition 
location for the high frequency case, however, led to discrepancies which may 
well be attributed to the pressure distribution obtained from inviscid-flow 
theory which does not consider any separation on the airfoil. This suggests 
the need for interaction between the viscous and inviscid flow calculations. 

Further investigations were conducted with the turbulence model of Cebeci 
and Smith which includes consideration of the transition region. It was found 
that the length of the transition region could be controlled by modifications 
to the model with consequent reduction of the transition region and of the 
separation bubble. Additional studies are required to quantify the effects, 
and, once again, should make use of interactive procedures. 
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Figure 8. Variation of transition location with angle of attack on the NAC A 
0012 airfoil at two reduced frequencies: (a) k = 0.01, (b) k = 
0.20, for case 2. The symbols denote experimental data of ref. 6 
and the solid lines the fairing done by the authors. 
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Figure 5. Variation of (a) displacenent thickness 6* and (b) local 

skjn-friction coefficient Cf for case ?., Eq. (19b), for k = 0.°0 
with transition corresponding to naxinun pressure peak. The arrow 
shows the upward stroke. 



(b) s/c 

Figure 6. Variation of (a) displacement thickness <$* and (h) local 

skin-friction coefficient Cf for case ?, Fq. (1%), for = D.?0 
with transition specified according to the experimental data, Fig 






rntire ft. The effect of the intermttency expr< ssior. on (a) displacement 
thickness <$*, and (h) local skin-friction coefficient r r -or 
case 2. Transition was calculated h w Michel's method. SnlH 1 
denote those computed with tie coefficient of 0 equal to l/l' >r,r ' 
dashed lines with 1/120. 
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